# This script performs link-predictive analysis of the PE network
# Last edited: 9/24/2013
# Edited by: Bruce Desmarais

# Read in network data
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/netList.RData")
# amat is the adjacency matrix recording the number of events both attend
# vdat is the vertex dataset containing the variables in DSMKdata

# functions to process vertex covariates
absdiffCov <- function(x){
	as.matrix(dist(cbind(x),diag=T,upper=T))
}

sameCov <- function(x){
	as.matrix(dist(cbind(x),diag=T,upper=T))==0
}


nodeCov <- function(x){
	matrix(x,length(x),length(x),byrow=T)+matrix(x,length(x),length(x),byrow=F)
}

overlapCov <- function(x,sep=":"){
	overMat <- matrix(0,length(x),length(x))
	for(i in 2:length(x)){
		for(j in 1:(i-1)){
			overMat[i,j] <- overMat[j,i] <- length(intersect(strsplit(x[i],split=sep)[[1]],strsplit(x[j],split=sep)[[1]]))
		}
	}
	overMat
}

listLengthCov <- function(x,sep=":"){
	xsp <- strsplit(x,split=sep)
	z <- unlist(lapply(xsp,length))
	nodeCov(z)
}

lowTri <- function(mat){
 	c(mat[lower.tri(mat)])
}

stDize <- function(mat){
	vecx <- c(mat)
	(mat - mean(vecx))/sd(vecx)
}


### Latent Space Models
require(latentnet)
estimLS <- list()
for(i in 1:10){
neti <- netList[[i]]$amat
vdati <- netList[[i]]$vdat

yi <- network(neti,dir=F)
el <- as.matrix(yi,matrix.type="edgelist")
set.edge.attribute(yi,"eval",neti[el])
nodeIdeo <- nodeCov(vdati$Dim1) # node-level effect of nominate
diffIdeo <- absdiffCov(vdati$Dim1) # absolute difference in nominate
sameParty <- sameCov(vdati$Party)   # same party
diffIdeoSP <- diffIdeo*sameParty # diffIdeo by same party
diffIdeoNSP <- diffIdeo*(1-sameParty) # diffIdeo by not same party
sameState <- sameCov(vdati$StateCode)   # same state
sameCohort <- sameCov(vdati$fow_seniority)   # same cohort
nodeMaj <- nodeCov(vdati$Party==median(vdati$Party))  # Dem node effect
nodeSeniority <- nodeCov(vdati$fow_seniority)  # seniority node effect
nodeElect <- nodeCov(vdati$Election==1)  # Election session effect
nodeCchr <- nodeCov(nchar(as.character(vdati$ChairT))>0) # Election session effect
commOverlap <- overlapCov(as.character(vdati$CommitteeT))   # Committee overlap
topicOverlap <- overlapCov(as.character(vdati$BillCodes))  # Bill topic overlap
nspon <- listLengthCov(as.character(vdati$BillCodes)) # Number of bills introduced by dyad
activity <- nodeCov(vdati$activity)  # solo press events


print(system.time(estimLS[[i]] <- ergmm(yi~euclidean(d=2)+latentcov(nodeIdeo)+latentcov(diffIdeoSP)+latentcov(diffIdeoNSP)+latentcov(sameParty)+latentcov(nodeMaj)+latentcov(nodeElect)+latentcov(nodeCchr)+latentcov(commOverlap)+latentcov(topicOverlap)+latentcov(nspon)+latentcov(activity)+latentcov(sameState)+latentcov(sameCohort)+latentcov(nodeSeniority), family="Poisson",response="eval")))

save("estimLS",file="~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/estimLS.RData")

print(i)

}

effNames <- c("edges","nodeIdeo","diffIdeoSP","diffIdeoNSP","sameParty","nodeMaj","nodeElect","nodeCchr","commOverlap","topicOverlap","nspon","activity","sameState","sameCohort","nodeSeniority")

## Ropeladder plots
coefList <- list()
for(i in 1:length(estimLS)){
	tabi <- summary(estimLS[[i]])$pmean$coef.table[,c(1,2,3)]
	coefList[[i]] <- tabi 
}


for(i in 1:length(effNames)){
filei <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/coef_",effNames[i],".pdf",sep="")
pdf(filei,height=3.05, width=5.75, family="Times", pointsize=12.25)
par(las=1,mar=c(3.5,4.5,.75,.5),cex.lab=2,cex.axis=1.75)
effMat <- NULL
for(j in 1:10){
	effMat <- rbind(effMat,coefList[[j]][i,])
	}	
	effMat <- as.matrix(effMat)
	plot(1:10, effMat[,1],type="n",ylab="",xlab="",xaxt="n",ylim=c(min(c(effMat)),max(c(effMat))))
	abline(h=0)
	abline(v=1:10,lty=2,col="grey60")
	segments(x0=1:10,y0=effMat[,2],y1=effMat[,3],lwd=3)
	points(1:10, effMat[,1],pch=19,cex=1.5)
	axis(1,lab=96:105,at=1:10,las=3)
	dev.off()
}


### Cosponsorship ###
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/netListC.RData")
library(latentnet)

estimLSc <- list()
for(i in 1:5){
neti <- netListC[[i]]$amat
vdati <- netListC[[i]]$vdat

yi <- network(neti,dir=F)
el <- as.matrix(yi,matrix.type="edgelist")
set.edge.attribute(yi,"eval",neti[el])
nodeIdeo <- stDize(nodeCov(vdati$Dim1)) # node-level effect of nominate
diffIdeo <- stDize(absdiffCov(vdati$Dim1)) # absolute difference in nominate
sameParty <- stDize(sameCov(vdati$Party))   # same party
diffIdeoSP <- stDize(diffIdeo*sameParty) # diffIdeo by same party
diffIdeoNSP <- stDize(diffIdeo*(1-sameParty)) # diffIdeo by not same party
sameState <- stDize(sameCov(vdati$StateCode))   # same state
sameCohort <- stDize(sameCov(vdati$fow_seniority))   # same cohort
nodeMaj <- stDize(nodeCov(vdati$Party==median(vdati$Party)))  # Dem node effect
nodeSeniority <- stDize(nodeCov(vdati$fow_seniority))  # seniority node effect
nodeElect <- stDize(nodeCov(vdati$Election==1))  # Election session effect
nodeCchr <- stDize(nodeCov(nchar(as.character(vdati$ChairT))>0)) # Election session effect
commOverlap <- stDize(overlapCov(as.character(vdati$CommitteeT)))   # Committee overlap
topicOverlap <- stDize(overlapCov(as.character(vdati$BillCodes)))  # Bill topic overlap
nspon <- stDize(listLengthCov(as.character(vdati$BillCodes))) # Number of bills introduced by dyad
activity <- stDize(nodeCov(vdati$activity))  # solo press events


print(system.time(estimLSc[[i]] <- ergmm(yi~euclidean(d=2)+latentcov(nodeIdeo)+latentcov(diffIdeoSP)+latentcov(diffIdeoNSP)+latentcov(sameParty)+latentcov(nodeMaj)+latentcov(nodeElect)+latentcov(nodeCchr)+latentcov(commOverlap)+latentcov(topicOverlap)+latentcov(nspon)+latentcov(activity)+latentcov(sameState)+latentcov(sameCohort)+latentcov(nodeSeniority),family="Poisson",response="eval",control=ergmm.control(group.deltas=.4,interval=500,sample.size=25000,burnin=1000000))))

save("estimLSc",file="~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/estimLSc1.RData")

print(i)

}

effNames <- c("edges","nodeIdeo","diffIdeoSP","diffIdeoNSP","sameParty","nodeMaj","nodeElect","nodeCchr","commOverlap","topicOverlap","nspon","activity","sameState","sameCohort","nodeSeniority")

estimLScCosp <- estimLSc

names(estimLScCosp) <- 96:105

## Ropeladder plots
coefList <- list()
for(i in 1:length(estimLScCosp)){
	tabi <- summary(estimLScCosp[[i]])$pmean$coef.table[,c(1,2,3)]
	coefList[[i]] <- tabi 
}


for(i in 1:length(effNames)){
filei <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/coef_",effNames[i],"C.pdf",sep="")
pdf(filei,height=3.05, width=5.75, family="Times", pointsize=12.25)
par(las=1,mar=c(3.5,4.5,.75,.5),cex.lab=2,cex.axis=1.75)
effMat <- NULL
for(j in 1:10){
	effMat <- rbind(effMat,coefList[[j]][i,])
	}	
	effMat <- as.matrix(effMat)
	plot(1:10, effMat[,1],type="n",ylab="",xlab="",xaxt="n",ylim=c(min(c(effMat)),max(c(effMat))))
	abline(h=0)
	abline(v=1:10,lty=2,col="grey60")
	segments(x0=1:10,y0=effMat[,2],y1=effMat[,3],lwd=3)
	points(1:10, effMat[,1],pch=19,cex=1.5)
	axis(1,lab=96:105,at=1:10,las=3)
	dev.off()
}

### Both Press Events and Cosponsorship Plots ###
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/estimLS.RData")
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/estimLSc.RData")

## Ropeladder plots

coefListC <- list()
for(i in 1:length(estimLSc)){
	tabi <- summary(estimLSc[[i]])$pmean$coef.table[,c(1,2,3)]
	coefListC[[i]] <- tabi 
}


coefList <- list()
for(i in 1:length(estimLS)){
	tabi <- summary(estimLS[[i]])$pmean$coef.table[,c(1,2,3)]
	coefList[[i]] <- tabi 
}


for(i in 1:length(effNames)){
filei <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/coef_",effNames[i],"Both.pdf",sep="")
pdf(filei,height=3.05, width=5.75, family="Times", pointsize=12.25)
par(las=1,mar=c(3.5,4.5,.75,.5),cex.lab=2,cex.axis=1.75)
effMat <- NULL
off <- .15
for(j in 1:10){
	effMat <- rbind(effMat,coefList[[j]][i,])
	}	
	effMat <- as.matrix(effMat)
effMatC <- NULL
for(j in 1:10){
	effMatC <- rbind(effMatC,coefListC[[j]][i,])
	}	
	effMatC <- as.matrix(effMatC)	
	plot(1:10, effMat[,1],type="n",ylab="",xlab="",xaxt="n",ylim=c(min(c(c(effMat),c(effMatC))),max(c(c(effMat),c(effMatC)))))
	abline(h=0)
	abline(v=1:10,lty=2,col="grey60")
	segments(x0=1:10-off,y0=effMat[,2],y1=effMat[,3],lwd=3)
	points(1:10-off, effMat[,1],pch=19,cex=1.5)
	segments(x0=1:10+off,y0=effMatC[,2],y1=effMatC[,3],lwd=3,col="grey55")
	points(1:10+off, effMatC[,1],pch=15,cex=1.5,col="grey55")
	axis(1,lab=96:105,at=1:10,las=3)
	dev.off()
}



